{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 第十章补充　Monte Carlo法：Monte Carlo积分\n",
    "\n",
    "作者：[王何宇](http://person.zju.edu.cn/wangheyu)\n",
    "\n",
    "[浙江大学数学科学学院](http://www.math.zju.edu.cn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function, division\n",
    "\n",
    "# 这句话是非标准的python，用于ipthon或jupyter这样的系统中，表示绘图即刻自动展开。\n",
    "%matplotlib inline\n",
    "\n",
    "# 这里把全部Warning过滤掉了. \n",
    "# 参见https://docs.python.org/2/library/warnings.html\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "from scipy.special import comb, perm\n",
    "from scipy import stats\n",
    "import numpy as np\n",
    "import sys\n",
    "import matplotlib.pyplot as plt\n",
    "np.random.seed(250)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "下面给一个一维MC积分的例子，\n",
    "$$\n",
    "I = \\int_a^b f(x) dx \\approx \\frac{b - a}{n} \\sum_{i = 1}^n f(X),\n",
    "$$\n",
    "其中$X \\sim U(a, b)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def obj_fun(x):\n",
    "    return np.sin(x)\n",
    "\n",
    "def MC_Int(a, b, times):\n",
    "    S = 0\n",
    "    X = np.random.rand(times)*(b - a)\n",
    "    for i in range(times):\n",
    "        S += obj_fun(X[i])\n",
    "    return (b - a) * S / times"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "MC_Int(0, 2*np.pi, 100000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对结构较为复杂的函数，更能体现优势。（注意在求面积时的例子有bug，已修改）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def obj_fun(x):\n",
    "    return 0.5 * np.sin(1/x) + 0.5\n",
    "\n",
    "MC_Int(0, 1, 100000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "为防止再次出bug，和经典方法做一个比较。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.integrate\n",
    "from numpy import sin\n",
    "f= lambda x:0.5*sin(1/x)+0.5\n",
    "i = scipy.integrate.quad(f, 0, 1)\n",
    "print (i)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "该结果不难推广到$m$维：\n",
    "$$\n",
    "\\int_{\\mathscr{J}^m} f(\\vec{x}) d\\vec{x} \\approx \\frac{\\lambda(\\mathscr{J}^m)}{n} \n",
    "\\sum_{i = 1}^n f(\\vec{X}) = \\frac{1}{n}\\sum_{i = 1}^n f(\\vec{X}),\n",
    "$$\n",
    "其中$\\vec{X}$是$\\mathscr{J}^m$上均匀分布的随机向量。下面是一个二维的例子："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def obj_fun2D(x, y):\n",
    "    return sin(x) * sin(y)\n",
    "\n",
    "def MC_Int2D(x0, x1, y0, y1, times):\n",
    "    S = 0\n",
    "    X = np.random.rand(times)*(x1 - x0)\n",
    "    Y = np.random.rand(times)*(y1 - y0)\n",
    "    \n",
    "    for i in range(times):\n",
    "        S += obj_fun2D(X[i], Y[i])\n",
    "    return (y1 - y0) * (x1 - x0) * S / times"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "MC_Int2D(0, 2*np.pi, 1, 3, 100000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "下面考虑方差估计，通过积分变换，将目标函数调整为$[0, 1]$上的简单函数："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def obj_fun(x):\n",
    "    return np.sin(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "按讲义第一个公式计算，我们需要在积分过程中保留全部函数计值，等待最终$\\bar{\\zeta}_n$计算出来之后，才能估计$\\mathrm{var} \\bar{\\zeta}_n$，这里为了和讲义一致，我们将随机数分布调整为$U(0, 1)$，而将积分变量调整到$[0, 1]$范围:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def MC_Int2(a, b, times):\n",
    "    S = 0\n",
    "    X = np.random.rand(times)\n",
    "    F = np.zeros(times)\n",
    "    for i in range(times):\n",
    "        F[i] = obj_fun(X[i] * (b - a))\n",
    "        S += F[i]\n",
    "    z = S / times\n",
    "    sigma2 = 0\n",
    "    for i in range(times):\n",
    "        sigma2 += (F[i] - z)**2\n",
    "    sigma2 = sigma2 / ((times - 1)*times)\n",
    "    return (z, sigma2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "MC_Int2(0, 2.0 * np.pi, 10000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们做一个检查，看有多少$\\bar{\\zeta}_n$落在$[\\zeta-\\frac{\\hat{\\sigma}_n}{\\sqrt{n}}, \\zeta+\\frac{\\hat{\\sigma}_n}{\\sqrt{n}}]$之间。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for T in range(100):\n",
    "    (z, v) = MC_Int2(0, 2.0 * np.pi, 1000)\n",
    "    e = np.sqrt(v)\n",
    "    z0 = z + e\n",
    "    z1 = z - e\n",
    "    plt.plot(T, z0, 'r.', T, z1, 'g.')\n",
    "plt.plot([0,100],[0,0],'--')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "如图，效果虽然不错（红点大多在虚线上而绿点大都在虚线下），但是这个估算要记录全部的$\\varphi(X)$的工作量过高。一个更常用的办法是："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def MC_Int(a, b, times):\n",
    "    S = 0\n",
    "    sigma2 = 0\n",
    "    X = np.random.rand(times)\n",
    "    for i in range(times):\n",
    "        F = obj_fun(X[i] * (b - a))\n",
    "        S += F\n",
    "        sigma2 += F**2\n",
    "    z = S / times\n",
    "    sigma2 = (sigma2 - times * z**2) / ((times - 1) * times)\n",
    "    return (z, sigma2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "MC_Int2(0, 2.0 * np.pi, 10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for T in range(100):\n",
    "    (z, v) = MC_Int2(0, 2.0 * np.pi, 1000)\n",
    "    e = np.sqrt(v)\n",
    "    z0 = z + e\n",
    "    z1 = z - e\n",
    "    plt.plot(T, z0, 'r.', T, z1, 'g.')\n",
    "plt.plot([0,100],[0,0],'--')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "效果差不多，但效率更高。我们之前分析过，直接用方差来做误差估计，会有一些缺陷（尽管在很多时候实际就这么用！），更准确的判断则需要做区间估计。这里回顾一下求面积时的区间估计。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def area_circle(times):\n",
    "    inside = 0\n",
    "    dots = np.random.rand(2, times)\n",
    "    for i in range(times):\n",
    "        x = dots[0, i] - 0.5\n",
    "        y = dots[1, i] - 0.5\n",
    "        if x * x + y * y < 0.25:\n",
    "            inside += 1\n",
    "    return inside / times"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "area_circle(100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "还是考虑一个圆面积的估算，我们抽样100次，得到一个估计值，那么区间估计要回答的问题是：真解在哪里？我们能否构建一个区间，覆盖住真解？完全覆盖不可能（总有人品爆炸的时候），那么能否给一个概率，称置信水平$\\delta$，也就是说能否给出一个区间，覆盖住真解的概率不小于$1 - \\delta$？"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "首先，这个问题同时也会涉及到参数$n$，也即采样次数。可以想像，在固定$\\delta$的前提下，如果$n$大，则区间可以小一些；如果$n$小，则区间会大一些。而在固定$n$的前提下，如果$\\delta$小，则区间范围会大；如果$\\delta$大，则区间范围会小。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "其次，我们考虑如何确定这个区间？作为一种统计估计，我们只能从统计量出发去估计。比如这里只有成功次数$S$，或者面积无偏估计$S/n$等信息（还有更高阶的）。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最后的问题就是我们如何理解统计的规律？如果它是服从概率论的某个定律，我们如何从中获取好处？"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "带着这些认识，我们重新复习一下区间估计这个关键问题。我们先将程序恢复成计数，也就是只输出成功次数$S$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def count(times):\n",
    "    inside = 0\n",
    "    dots = np.random.rand(2, times)\n",
    "    for i in range(times):\n",
    "        x = dots[0, i] - 0.5\n",
    "        y = dots[1, i] - 0.5\n",
    "        if x * x + y * y < 0.25:\n",
    "            inside += 1\n",
    "    return inside"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "count(20) # 多抽几次，让它等于14"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上面是做一次$n = 20$的估计得到的成功次数，现在是直面人品的时候了！在这一次$n = 20$的试验中，我们最有可能得到的结果是多少？这个问题的回答并不唯一，因为我们拥有的信息可能不一样。有些时候，我们对整个抽取过程了解的清楚一点，比如我们知道具体的抽样过程；有些时候，可能只知道这是个独立的抽取过程；甚至有时连独立性也要先判定。所以一切概率分析的出发点，是要先做个假设！现在我们找个好日子过，我们清楚的知道，这是对一个均匀分布的随机投点，是否会落入一个待测面积的区域的判定。这样，由于待测区域是固定不变的，落入期间的概率也是一个定值：\n",
    "$$\n",
    "p = \\frac{1}{\\lambda(\\mathscr{R})}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "当然这个$p$我们并不知道（这就是我们要估计的目标！），但我们清楚，这整个抽取过程就是一个$n = 20$的二项分布，也即$S$服从概率分布：\n",
    "$$\n",
    "p_i = P\\{S = i\\} = \\binom{n}{i} p^i(1 - p)^{n - i}.\n",
    "$$\n",
    "相应的累积分布为：\n",
    "$$\n",
    "F(i) = \\sum_{j = 0}^i\\binom{n}{j} p^j(1 - p)^{n - j}.\n",
    "$$\n",
    "\n",
    "可以画出来看看："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.special import comb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def binomial_CDF(n, mu, i):\n",
    "    F = 0\n",
    "    for j in range(i + 1):\n",
    "        F += comb(n, j) * mu**j * (1 - mu)**(n - j)\n",
    "    return F"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 20\n",
    "S = np.zeros(n + 1)\n",
    "F1 = np.zeros(n + 1)\n",
    "F2 = np.zeros(n + 1)\n",
    "F3 = np.zeros(n + 1)\n",
    "for i in range(n + 1):\n",
    "    S[i] = i\n",
    "    F1[i] = binomial_CDF(n, 0.1, i)\n",
    "    F2[i] = binomial_CDF(n, 0.5, i)\n",
    "    F3[i] = binomial_CDF(n, 0.9, i)\n",
    "plt.bar(S, F1)\n",
    "plt.bar(S, F2)\n",
    "plt.bar(S, F3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "现在我们假设$S$服从二项分布，在一组试验中，得到了$S = 16$。我们来估计$p$的可能存在区间。先随便取一个$p=0.7$："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 20\n",
    "S = np.zeros(n + 1)\n",
    "F1 = np.zeros(n + 1)\n",
    "for i in range(n + 1):\n",
    "    S[i] = i\n",
    "    F1[i] = binomial_CDF(n, 0.7, i)\n",
    "plt.plot(S, F1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们看到，$S$最有可能大概是14左右，这个其实是确定，已知的，因为我们现在明确知道$n = 20$，$p = 0.7$，于是最大可能就是$S = 14$。但这里$p = 0.7$是我们乱猜的。我们事实上不可能知道真正的$p$。这里我们说的是，如果$p$是$0.7$，那么最有可能的是$14$。然而我们现在得到的是$16$，是不是就错了？当然不是。事实上在$p = 0.7$前提下，$S$恰好等于$14$，只是相对最有可能，绝对数值仍然不是很大："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "binomial_CDF(20,0.7,14) - binomial_CDF(20,0.7,13)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "只有不到两成机会。现在来换一个角度看，$S$有多大的机会严格小于$16$?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "binomial_CDF(20,0.7,15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "看来$p = 0.7$似乎是个很糟糕的估计，因为在这个估计下，$S < 16$的可能性不是个小量。但问题的关键来了，**我们究竟是希望这个数再大点好？还是再小点好？**我们现在估计的可不是$S$，我们在估计$p$，我们想得到$p$的**最有可能**范围，因此，我们首先应该去掉$p$的最不可能范围。对于一个**已知**的统计结果$S = 16$，我们现在看到实际上如果$p = 0.7$，则$S < 16$的可能性还有$0.76$，这个不够大，也就是说$p = 0.7$还是很有可能发生的，因为在这个前提下，会有大于$0.2$的概率，出现$S \\geq 16$的情况。而我们希望，首先确定一个小概率，比如$\\frac{\\delta}{2} = 0.05$，使得发生$S \\geq 16$的可能性小于$\\frac{\\delta}{2}$。那么我们怎么调整$p$？往大调还是往小调？观察图像，显然$p$应该再小点。也即我们寻找$p$，使得\n",
    "$$\n",
    "\\alpha_1 = 1 - F_{15}(n, p) = \\frac{\\delta}{2},\n",
    "$$\n",
    "这样，当$p > \\alpha_1$时，出现$S \\geq 16$的概率就会大于$\\alpha_1$，换言之，**去掉的部分，发生的可能性小于$\\alpha_1$**。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a1 = 1 - binomial_CDF(20,0.7,15)\n",
    "print(a1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这个$p = 0.7$显然不对，我们需要解上面这个方程，二分法考虑一下："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def func(p, S, n, alpha):\n",
    "    return binomial_CDF(n, p, S) - alpha\n",
    "\n",
    "def bisection_root(alpha, S, n, p1, p2, err, tol):\n",
    "    while (True):\n",
    "        f1 = func(p1, S, n, alpha)\n",
    "        if np.abs(f1) < tol:\n",
    "            return p1\n",
    "        p = (p1 + p2) * 0.5\n",
    "        fz = func(p, S, n, alpha)\n",
    "        if np.abs(fz) < tol:\n",
    "            return p\n",
    "        if fz * f1 < 0:\n",
    "            p2 = p\n",
    "        else:\n",
    "            p1 = p\n",
    "        p = (p1 + p2) * 0.5\n",
    "        if (np.abs(p1 - p2) < err):\n",
    "            return p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "EPS = 1e-7 #精度\n",
    "p1 = bisection_root(0.995 + EPS, 15, 20, 0, 1, EPS, EPS)\n",
    "print(p1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "注意上面加一个误差界以确保$p < p_1$时$S \\geq 16$的概率能严格小于$\\frac{\\delta}{2}$。检查一下："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "1 - binomial_CDF(20,p1,15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$p$稍微大一点，$S >= 16$发生的概率就超过$\\frac{\\delta}{2}$了:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "1 - binomial_CDF(20, p1 + EPS * 2, 15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "换言之，既然现在我们观察到$S = 16$，又要求有$1 - \\frac{\\delta}{2}$的概率正确，那么，$p$必须大于$p_1$。由此我们找到了区间估计的第一个界。对称地，我们去寻找第二个界。$S > 16$的可能性要足够小！"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p2 = bisection_root(0.005 - EPS, 16, 20, 0, 1, EPS, EPS)\n",
    "print(p2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "检查一下，如果$p > p2$，那么$S \\leq 16$的概率就会严格小于："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "binomial_CDF(20, p2, 16)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "不能再稍微小一点："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "binomial_CDF(20, p2 - EPS * 2, 16)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "综合起来，对$S = 16$这个观测结果而言，$p$落在$(p_1, p_2)$的可能性严格大于$1 - \\delta$。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "总结一下，实际上区间估计，不是估计！它是从一个实际统计量出发，在先验假设（服从二项分布）的前提下，直接拎出两个分布函数，一个是$F_{S - 1}(n, p)$，一个是$F_S(n, p)$，然后在偏左的$F_{S - 1}$中，直接砍掉$[0, \\alpha_1]$部分，一刀砍在$p_1$；而在偏右的$F_S$中，砍掉$[\\alpha_2, 1]$部分，一刀砍在$p_2$，从而确保在$p \\in (p_1, p_2)$的可能性严格大于$\\alpha_2 - \\alpha_1$。一般我们取$\\alpha_1 = \\frac{\\delta}{2}$，$\\alpha_2 = 1 - \\frac{\\delta}{2}$，于是$P\\{p_1 < p < p_2\\} > 1 -\\delta$。显然，从不同的假设出发，我们的先验分布未必是二项分布，也可以是正态分布（$n$充分大），或其他合理的假设。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
